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ABSTRACT 


The  mechanical  and  rheological  properties  of  an  MR  fluid  depend  on  the  induced 
microstructure  of  the  imbedded  ferrous  particles.  When  subject  to  an  externally  applied 
magnetic  field  these  particles  magnetize  and  align  themselves  in  chains  parallel  to  the 
applied  field.  The  microstructure  of  these  chains  is  a  function  of  several  parameters 
including  particle  size,  applied  magnetic  field  strength,  and  viscosity  and  velocity  of  the 
surrounding  fluid.  This  thesis  will  create  a  model  from  a  first  principle  approach  to 
accurately  predict  the  micro  structure  in  a  variety  of  different  situations.  The  model 
investigated  assumes  the  particles  become  magnetic  dipoles  upon  the  application  of  the 
magnetic  field  and  that  particle  interaction  is  due  solely  to  dipole-dipole  interaction.  Due 
to  the  inherently  small  size  of  the  particles,  drag  is  modeled  using  Stokes’  drag.  This 
mathematical  model  will  be  used  to  create  a  computer  simulation  to  visualize  and  analyze 
the  subsequent  transient  micro  structures  formed.  The  model  will  assume  a  constant 
magnetic  field  applied  (i.e.,  no  spatial  or  time  gradients)  and  that  the  effects  of  this  field 
are  felt  instantaneously. 
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I.  INTRODUCTION 


A.  CONTROLLABLE  FLUIDS 


Magnetorheological  (MR)  and  electrorheological  (ER)  fluids  are  a  class  of 
“smart”  materials  that  are  characterized  by  their  ability  to  reversibly  transform  from 
liquid  state  to  a  Bingham  solid.  They  are  fluids  that  have  either  magnetically  permeable 
(or  electrically  conductive)  microscopic  particles  suspended  in  them.  The  transformation 
from  liquid  state  to  Bingham  solid  occurs  by  the  application  of  a  magnetic  (or  electric) 
field  to  the  fluid.  This  magnetic  (or  electric)  field  causes  the  suspended  particles  to  align 
in  chains  along  the  field  lines  in  a  manner  to  reduce  the  overall  energy  of  the  field.  The 
existence  of  these  chains  changes  many  bulk  properties  of  the  fluid.  Of  practical  interest 
is  the  change  in  viscosity  which  causes  the  fluid  to  behave  like  a  Bingham  solid. 

The  Bingham  model  used  for  modeling  MR  fluids  relates  the  total  shear  stress  z 
to  the  shear  rate  y  and  H  (magnitude  of  applied  magnetic  field)  according  to  the  equation 


T  = 


Ty(H)  +  1) 


7 


sgn  (y) , 


where  rv  (//)  is  a  yield  stress  that  is  a  function  of  the  applied  magnetic  field  and  i)  is  the 

composite  bulk  viscosity  of  the  fluid  [1].  This  equation  is  phenomenological  in  nature 
where  the  values  in  the  equation  are  determined  experimentally  instead  of  being  deduced 
from  a  first  principle  approach. 

Recently  it  has  been  found  that  a  more  detailed  approach  to  predicting  the 
behavior  of  MR  fluids  has  become  necessary  due  to  the  limitations  of  the  above  approach 
[1].  First,  the  Bingham  model  is  a  macro  scale  approach  (the  fluid  and  particles  are 
treated  as  a  single  continuum  instead  of  a  composite  system)  with  no  differentiation  with 
particle  level.  The  coupling  between  mechanical  behavior  and  the  magnetic  field  takes 
place  at  the  particle  level  and  is  governed  by  first  principles  (conservation  of  momentum, 
Maxwell’s  Laws,  etc.).  The  Bingham  model  then  is  limited  to  a  narrow  range  of 
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applicability  commensurate  with  the  experimental  data.  Second,  the  Bingham  model 
tends  to  be  inaccurate  at  low  value  of  stress.  In  current  applications  where  MR  devices 
are  used  as  feedback  controls,  low  value  stresses  are  important  and  the  above  model 
proves  unsatisfactory.  Third,  the  Bingham  model  is  only  applicable  to  1-D  simple  shear 
flows  with  a  transverse  magnetic  field  applied.  This  is  inadequate  for  multi-degree  of 
freedom  MR  devices  that  are  currently  being  designed.  These  reasons  encourage  a 
different  model  to  be  developed  that  is  based  on  first  principles. 

B.  INDUSTRIAL  APPLICATIONS 

The  American  inventor,  Willis  Winslow,  was  the  first  to  recognize  how  to  create 
a  smart  fluid  using  these  principles  [2,  3],  He  did  much  of  the  initial  pioneering  work  on 
ER  fluids  in  the  1940s.  Later,  Jacob  Rabinow  investigated  the  same  phenomenon  using  a 
magnetic  field  for  use  in  a  magnetic  field  clutch  [2]  and  is  considered  the  first  to  develop 
MR  fluids.  Although  their  works  were  conducted  over  half  a  century  ago,  it  has  only 
been  recently  that  the  use  of  these  smart  fluids  has  become  more  common  in  industrial 
applications.  This  is  due  primarily  to  the  stability  and  durability  requirements  of  modern 
designs  [3]. 

Today  the  uses  for  these  smart,  controllable  fluids  are  numerous  and  varied.  One 
primary  use  is  in  hydraulic  dampers  and  brakes.  Because  of  the  ability  to  rapidly  change 
the  working  fluid  viscosity,  one  has  the  ability  to  change  the  damping  coefficient  in 
dampers  to  give  much  better  dynamic  response  and  control.  Other  applications  include 
better  feedback  to  control  items  such  as  joysticks,  responsive  personnel  armor,  and  MR 
polishing  machines  [4], 
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c. 


TYPICAL  MR  FLUID  ARRANGEMENT 


A  typical  arrangement  for  an  industrial  application  of  a  MR  fluid  is  shown  below. 


MR  Fluid 

Electromagnet 

* 

v  A 

r  i 

f  ' 

f  ' 

f  ^ 

f 

Magneticr 

Field  Lines 

/s'  Electromagnet 

Figure  1 .  MR  Fluid  Device  Arrangement 

The  MR  fluid,  consisting  of  a  carrier  fluid  (usually  a  silicone  oil)  and  the 
suspended  particles  (typically  fine  ferrous  particles),  is  sandwiched  in  a  small  gap 
between  two  electromagnets.  These  magnets,  when  energized,  create  a  magnetic  field 
perpendicular  to  the  flow  of  the  MR  fluid  which  causes  the  imbedded  particles  to  form 
chains  parallel  to  the  applied  field.  The  dynamic  response  of  these  particles  in  both  a 
static  fluid  and  a  moving  fluid  is  investigated  in  this  thesis. 
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II.  MAGNETIC  FORCE 


A.  DIPOLE  MODEL 

When  the  magnetic  field  is  applied  to  the  MR  fluid,  the  ferromagnetic  particles 
become  magnetized  and  interact  with  the  field  and  with  each  other.  The  exact  solution  to 
these  interactions  is  difficult  (if  not  impossible)  and  involves  the  integration  of  Maxwell’s 
stress  tensors  across  the  entire  volume  of  the  magnetized  solution.  Instead,  some 
simplifying  assumptions  need  to  be  made  in  order  to  allow  an  easily  calculable  analytical 
solution  without  sacrificing  accuracy. 

The  first  assumption  to  be  made  is  to  use  a  dipole  model  for  the  magnetized 
particles.  This  assumes  that  the  particles  are  dumbbell  in  shape,  with  length  L.  One  end 
contains  the  positive  (North)  pole  and  the  other  contains  the  negative  (South)  pole  of  the 
magnet  as  shown  in  Figure  2.  The  magnitude  of  the  pole  strength  is  denoted  by  q  and 
arises  because  the  applied  magnetic  field  magnetizes  the  particle. 


L 


Figure  2.  Dumbbell  Shape 

When  the  above  dumbbell  is  placed  in  a  magnetic  field  H,  it  experiences  a  torque 
about  its  center  as  described  by  the  formula  r  =  LqH  sin(#) ,  where  6  is  the  angle  between 
the  magnetic  field  and  dumbbell.  The  quantity  L*q  is  given  a  special  name,  the  magnetic 
moment,  and  is  denoted  by  m.  The  model  used  in  this  thesis  uses  the  dipole  model  and  is 
defined  by  determining  the  value  of  m  in  the  limit  where  L  goes  to  zero  but  the  torque 
remains  finite.  This  is  the  case  of  magnetic  spheres  which  are  used  in  most  MR 
applications. 
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The  magnitude  of  the  dipole  moment  determines  the  interaction  of  the  particles 
with  the  external  magnetic  field  and  the  interaction  of  the  particles  with  each  other.  It  is  a 
function  of  the  magnitude  of  the  applied  external  field  (H  with  units  Amp/m),  the  volume 
of  the  particles,  and  the  magnetic  permeability  of  both  the  particles  and  the  surrounding 
fluid  [5].  Specifically  it  is  given  by  the  relation 

r  \ 


m  =  (4  nas') 


Mp-Mf 


H , 


loc 


(1) 


where  a  =  radius  of  the  particle  (meters),  jup  =  magnetic  permeability  of  the  particle 
(henry/meter),  /Uf  =  magnetic  penneability  of  the  fluid,  and  Hioc  =  magnetic  field  at  dipole 
location  (amp/meter). 

Several  assumptions  need  to  be  made  when  using  the  above  formula.  The 
presence  of  a  dipole  alters  the  magnetic  field  in  its  vicinity  (this  is  what  causes  particles 
to  interact  with  each  other)  and  this  implies  that  the  value  of  Hioc  needs  to  be  calculated  at 
every  point.  However,  this  variation  is  assumed  to  be  negligible  when  calculating  the 
magnetic  dipole  since  the  external  fields  applied  are  relatively  large  (on  the  order  of  200 
kA/m)  and  the  variations  caused  by  the  dipoles  are  several  orders  of  magnitude  smaller. 
Therefore  Hioc  is  assumed  to  be  equal  to  the  applied  magnetic  field.  The  second 
assumption  concerns  the  value  of  iup.  Because  the  particles  are  ferrous,  /up  is  not  a 
constant  but  varies  with  the  applied  magnetic  field.  However  since  the  range  of  the 
applied  field  is  small,  often  a  fixed  value,  an  average  value  of  /up  is  used  based  on  the 
values  of  the  applied  field. 


B.  FORCE  ON  A  DIPOLE 


The  force  on  a  dipole  in  a  magnetic  field  is  given  by  the  product  of  the  dipole 
moment  and  the  gradient  of  the  magnetic  field  as  given  by  the  expression 


F 


dH 

=  m - 

dr 


(2) 
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where  Fr  is  the  force  in  some  arbitrary  direction  r  and  - —  is  the  spatial  derivative  in  the 

dr 


r  direction.  One  method  that  presents  itself  in  determining  the  forces  is  to  simply 
calculate  the  magnetic  field  at  every  point.  Theoretically  this  could  be  done  by 
calculating  the  external  applied  field  and  modifying  it  by  the  perturbations  caused  by  the 
presence  of  the  dipoles.  Since  the  location  of  the  dipoles  constantly  changes,  this 
calculation  would  have  to  be  perfonned  at  every  time  step.  Using  this  calculation  (which 
would  have  to  be  perfonned  numerically)  the  gradient  at  every  point  could  be  calculated 
and  then  the  force  on  every  particle  could  be  determined.  In  reality  this  calculation  is 
difficult  to  perform,  requires  specific  algorithms  for  determining  the  field  and  the 
gradients,  and  requires  massive  computing  power.  A  more  simplistic  approach  was  used 
for  this  model. 

The  first  assumption  for  a  more  simplistic  approach  is  that  the  applied  magnetic 
field  is  uniform.  This  assumption  is  valid  since  the  applied  magnetic  field  is  enacted 
rapidly  (assumed  instantaneous),  does  not  vary  with  time,  and  the  fluid  gap  is  very  small 
compared  to  the  surface  area  over  which  the  field  is  enacted.  Since  the  magnetic  field  is 
assumed  uniform  in  space  and  time  there  is  no  gradient  and  the  particles  experience  no 
net  force  due  to  the  external  field.  The  only  magnetic  force  the  particles  experience  is 
due  to  their  mutual  interactions. 

Consider  two  magnetic  dipoles  of  identical  strength  at  arbitrary  positions  r.  and 
r. .  A  magnetic  field  H0  is  applied  parallel  to  the  z  axis.  The  force  between  the  two 
dipoles  is  given  by 


r  —r. 

<  j 


0r  is  the  angle  from  the  z 


where  ftj  is  the  force  on  particle  i  from  particle  j,  r..  = 
axis  and  r  ,  er  is  a  unit  vector  parallel  to  r..  ,  and  eg  is  a  unit  vector  parallel  to 
er  x (e,  xHq)  [6].  This  is  shown  in  the  below  figure. 
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Figure  3.  Relationship  Between  Two  Particles 


This  equation  models  the  interaction  between  dipoles  and  it  is  useful  to  examine 
this  equation  quantitatively  to  obtain  a  feel  for  the  dynamics  of  the  particles.  By 
combining  equations  (1)  and  (2),  it  becomes  apparent  that  the  interaction  force  is 
proportional  to  the  square  of  the  applied  field,  proportional  to  the  square  of  the  particles 
volume,  and  the  direction  of  the  force  is  a  function  of  the  relative  location  of  the  two 
particles.  This  last  item  is  what  causes  the  particles  to  form  stable  chains  when  the 
magnetic  field  is  applied.  Examine  only  the  radial  term  in  equation  (3).  If  0tj  is  less  than 

-54.6  degrees,  the  particles  tend  to  attract.  Otherwise  they  tend  to  repel. 

It  is  more  useful  to  transform  equation  (3)  into  Cartesian  coordinates.  Using  the 
same  x,  y,  and  z  directions  as  shown  in  Figure  (2)  and  defining  Q  =  3 m2jur ,  the  x,  y,  and  z 
components  are  given  as 
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ru 
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where /y*  is  the  x  component  of fj,fjy  is  the  y  component  of fj,  and/yz  is  the  z  component 
of fij.  The  derivations  for  the  above  equations  are  attached  in  Appendix  A. 

In  a  suspension  of  N  particles  each  with  an  assumed  induced  magnetic 
dipole  moment  of  m,  the  total  magnetic  force  due  to  dipole  interaction  on  a  particle  i  is 
the  sum  of  the  contributions  of  all  of  the  other  particles  in  the  suspension.  In  algebraic 
form 


(7) 

(8) 

(9) 


where  FMix  is  the  total  magnetic  force  on  particle  i  in  the  x  direction,  FUiy  is  the  total 

magnetic  force  on  particle  i  in  the  y  direction,  and  Fm_  is  the  total  magnetic  force  on 

particle  i  in  the  z  direction.  To  determine  the  dynamic  behavior  of  the  particles  in  the 
fluid,  these  equation  are  calculated  at  every  time  step,  the  deviation  in  the  current  position 
of  the  particles  are  calculated,  the  values  of  the  forces  are  recomputed  at  the  next  time 
step  based  on  the  particles’  new  position,  and  the  process  is  repeated  until  the  end  of  the 
computational  time. 
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III.  OTHER  FORCES 


A.  DRAG  FORCE 

The  drag  force  on  a  spherical  particle  moving  in  a  viscous  fluid  is  a  function  of 
the  pressure  difference  across  the  sphere  (fonn  drag)  and  the  surface  shear  stress  (viscous 
drag).  In  general  this  expression  can  be  complicated  to  solve.  In  the  specific  case  of  the 
small  particles  used  in  MR  fluids  a  number  of  simplifying  assumptions  can  be  made  to 
more  easily  determine  this  drag  force.  The  flow  can  be  assumed  to  be  laminar  due  to  the 
small  clearances  between  the  electromagnets  which  the  fluid  flows  between  and  the  small 
velocities  analyzed  in  this  thesis.  Another  simplifying  case  arises  due  to  the  small 
particle  size  (-micro  meter)  which  implies  that  the  Reynolds  number  based  on  diameter 
(Rcd)  is  less  than  1.  Both  of  these  assumptions  allow  the  viscous  drag  force  to  be 
modeled  by  the  well  understood  Stokes’  drag  which  is  given  by 

dr 

Fr=bmja —  (10) 

dt 

where  Fr  is  the  drag  force  in  the  r  direction,  77  is  the  viscosity  of  the  fluid  and  a  is  the 
radius  of  the  particle. 

B.  GRAVITY 

The  force  of  gravity  is  neglected  in  this  model  based  on  the  fact  that  the 
gravitational  force  that  would  tend  to  make  the  particles  settle  is  a  much  weaker  force 
than  the  magnetic  force  that  acts  between  the  dipoles.  This  is  obviously  not  true  when  no 
magnetic  field  is  applied,  but  the  gravitational  settling  is  ignored  by  assuming  that  the 
suspension  is  thoroughly  mixed  before  the  application  of  the  field  and  that  the  particles 
are  randomly  distributed  in  the  carrier  fluid. 
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c. 


BROWNIAN  MOTION 


Brownian  motion  is  characterized  by  the  random  walk  of  particles  in  a  fluid  due 
to  the  bombardment  of  molecules.  In  detennining  if  Brownian  motion  should  be 
considered  in  any  model  of  MR  fluids,  this  effect  should  be  compared  to  other  effects 
which  determine  the  dynamic  behavior  of  the  particles.  When  there  is  no  applied 
magnetic  field  this  is  not  the  case,  however  by  assuming  that  the  fluid  is  mixed 
immediately  before  the  field  is  applied  would  negate  any  effects  of  Brownian  motion. 
Once  the  field  is  applied  the  relative  effect  of  Brownian  motion  compared  with  the 
magnetic  forces  can  be  detennined  by  analyzing  the  coupling  constant  X  which  is  defined 
as  the  ratio  of  the  interaction  energy  of  two  dipoles  in  contact  and  the  thermal  energy  [6]. 
Specifically 

z=Kl  (11) 

KT  9i„T 

where  /u0  is  the  magnetic  permeability  in  a  vacuum,  a  is  the  particle  radius,  II  is  the 
magnetic  field,  x  is  the  magnetic  susceptibility  of  the  particle,  kb  is  the  Boltzmann 

constant  and  T  is  the  absolute  temperature.  In  all  cases  modeled  in  this  thesis  X  »  1  and 
Brownian  motion  is  ignored. 

D.  REPULSIVE  FORCES 


The  particles  themselves  and  any  walls  that  physically  constrain  the  MR  fluid  are 
modeled  as  hard  surfaces.  Therefore  a  fictitious  repulsive  force  must  be  modeled  to 
ensure  that  a  particle  in  physical  contact  with  either  a  wall  or  another  particle  behaves  as 
hard.  The  characteristics  of  this  force  are  such  that  when  two  particles  touch  (the 
distance  between  two  dipoles  is  2 a  apart)  the  repulsive  force  exactly  balances  the 
attractive  force  between  the  dipoles,  and  when  the  distance  between  the  dipoles  is  greater 
than  2 a  the  force  is  negligibly  small.  The  proposed  form  of  this  force  is  given  below 


k2 

f  ..=K,e 

J  rep  ,ij  1 


(12) 
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where  K\  and  Kj  are  constants  to  be  determined.  The  exponential  term  was  chosen  to 
give  a  function  that  rapidly  decays  as  r(/-  increases  and  is  a  commonly  used  mathematical 
model  for  these  types  of  interactions  [7]. 

To  determine  K\,  apply  the  condition  that  when  r-j  is  equal  to  2 a  then  the  repulsive 
force  must  equal  the  attractive  force  between  the  dipoles.  From  equation  (3)  the  dipoles 
attract  when  6f  is  less  than  -54.6  degrees  and  the  attractive  force  also  causes  0..  to  tend  to 

zero.  This  is  what  causes  the  particles  to  align  in  chains  that  characterize  the  MR  fluid. 
Assume  that  the  particles  will  touch  when  6{j  is  small.  From  Figure  2,  this  implies  that 

Xj  and  y..  are  negligibly  small.  Appling  this  condition  to  equations  (4-6)  gives 

-4  =  °’  -4  = 0 > and  4=4 

V 


3- 


A 


When  the  particles  are  touching  z..  =  r  =  2 a  which,  when  combined  with  the 
above  equation,  gives 


2Q 

(2 A  ' 


Combining  this  result  with  equation  (9)  gives  a  value  for  Kl  = 


Q 


8a 


4  ' 


The  constant  Ki  is  detennined  by  a  much  less  rigorous  means.  It  must  be  negative 
to  give  a  decaying  characteristic  and  its  magnitude  is  selected  by  a  trial  and  error 
approach.  On  one  hand  a  high  magnitude  gives  a  steeper  decay  which  is  advantageous 
since  this  more  closely  approximates  reality.  However,  if  the  value  is  too  large,  the 
repulsive  term  can  become  extremely  large  for  small  distances  and  leads  to  numerical 
instabilities.  A  value  of  K2  =  -12  was  chosen  as  a  balance  between  these  two  competing 

factors  based  on  numerical  experiments.  This  gives  a  steep  decay  while  allowing  a  more 
manageable  time  step. 
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Using  the  values  of  K\  and  /C  and  transforming  equation  (9)  into  Cartesian 
coordinates  gives  the  following  expressions  for  the  repulsive  force  on  a  particle  i  from 
particle  j  in  the  x,  y,  and  z  directions  as 


(13) 


(14) 


(15) 


For  the  physical  interaction  with  the  walls  of  the  container,  a  similar  approach 
was  taken,  and  the  equation  developed  for  the  interaction  of  a  particle  with  the 
floor/ceiling  is  given  as 


(16) 


where  fWiz  is  the  force  on  particle  i  from  the  wall  in  the  z  direction,  z.  is  the  absolute 
distance  from  the  bottom  boundary  to  particle  i  in  the  z  direction,  and  H  is  the  total  height 
of  the  volume  (not  to  be  confused  with  the  use  of  H  elsewhere  as  the  magnetic  field 
strength).  Equations  identical  to  (16)  are  used  for  the  horizontal  boundaries  with  the 
substitutions  for  the  particles  x  and  y  positions  (instead  of  z. )  and  the  length  and  width  of 
the  contaimnent  area  (instead  of  H). 

E.  INTERACTION  WITH  ELECTROMAGNET 

Up  to  this  point  the  discussion  of  the  physics  of  the  interactions  of  the  particles  in 
a  MR  fluid  is  exactly  the  same  as  if  it  was  an  ER  fluid  (replace  the  electromagnet  with 
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charged  parallel  plate  conductors  and  replace  some  of  the  magnetic  constants  with  their 
electrical  equivalents).  A  primary  difference  between  the  two  arises  in  the  physics  of  the 
interaction  between  the  electromagnet  (MR)  and  the  interaction  with  the  charged 
conducting  plate  (ER).  In  the  latter  case,  the  charged  plates  induce  an  electric  dipole 
(exactly  analogous  to  the  electromagnets  inducing  a  magnetic  dipole),  but  the  electric 
dipoles  interact  with  the  plates  in  an  easily  understood  manner.  The  presence  of  the 
electric  dipoles  themselves  will  induce  a  current  distribution  on  the  plates  and  then  these 
dipoles  are  electrically  attracted  to  the  plates  because  of  this  current  distribution.  A  well 
documented  manner  to  calculate  this  interaction  is  by  the  method  of  image  charge  [5]. 
Basically  the  interaction  of  a  dipole  that  is  a  distance  L  from  the  plate  is  identical  to 
assuming  there  is  an  infinite  number  of  equal  dipoles  on  the  other  side  of  the  plate  at 
distances  nL  where  n=l,2,3....  Therefore  an  electric  dipole  will  interact  with  the 
conducting  plate,  specifically  will  be  attracted  to  the  plate  and  attach  itself  to  the  plate.  If 
there  are  multiple  dipoles,  they  will  form  chains,  and  the  chains  will  anchor  themselves  to 
the  plate  and  behave  as  if  the  chain  was  infinitely  long.  This  is  what  allows  an  ER  to 
have  a  shear  stress;  the  chains  are  anchored  to  the  plate. 

There  is  no  analogy  in  the  MR  case.  There  is  no  such  thing  as  a  magnetic  current 
produced  at  the  boundary  of  the  electromagnet  that  would  allow  the  use  of  the  dipole 
image  method  to  detennine  the  interaction  of  the  chain  with  the  magnet  [6].  Another  way 
to  look  at  this  is  to  consider  a  single  magnetic  dipole  between  the  magnets.  Assuming  a 
constant  magnetic  field,  the  dipole  would  not  be  attracted  to  either  of  the  magnets.  There 
seems  to  be  nothing  to  lock  the  chain  in  place  and  therefore  an  MR  fluid  would  not  be 
able  to  have  a  shear  stress.  Experimentally,  this  is  not  true.  The  chains  do  become 
locked. 

There  are  two  reasonable  theories  as  to  how  this  locking  occurs.  The  first  is  to 
question  the  assumption  that  the  field  is  uniform.  Away  from  the  edges  of  the  magnets, 
due  to  the  small  gap  between  the  magnets,  it  is  safe  to  assume  a  uniform  field.  When  the 
magnets  are  close  together  the  field  lines  away  from  the  edges  do  not  spread  out  and 
consequently  there  is  no  gradient.  However,  at  the  edges  of  the  magnets,  this  is  not  true. 
The  fields  bulge  outward  and  tend  to  wrap  around,  causing  large  gradients.  One  proposal 

15 


is  that  away  from  the  edges  of  the  magnets,  the  chains  are  free  to  move  (are  not  locked  to  the 
magnet),  but  as  the  bulk  fluid  flow  sweeps  them  toward  the  edge  of  the  magnet,  they  become 
locked  in  this  area  of  high  field  gradients  and  effectively  form  a  lattice  type  wall.  Other 
chains  being  swept  along  will  then  build  up  behind  this  lattice  wall. 

Another  theory  approach  is  to  again  to  question  the  uniformity  of  the  field,  this  time 
at  the  fluid/magnet  interface.  To  explain  this  effect  requires  the  use  of  Maxwell’s  laws.  The 
equation  of  specific  use  is 

V»B  =  0  (17) 

where  B  is  magnetic  flux  density  (Tesla).  The  relationship  between  H  and  B  is  given  by 

B  =  /jH  (18) 

where  /u  is  the  magnetic  permeability  of  the  substance  through  which  H  exists.  Using 
equations  (17)  to  solve  for  the  normal  component  of  B  across  the  discontinuity  between  the 
magnet  and  the  fluid  (there  are  no  tangential  components)  implies 

nl»(/i2H2  -  H] )  =  0  or  jU2H2  =  //,//, 

where  the  subscripts  refer  to  the  magnetic  permeability  and  magnetic  field  of  the  magnet  and 
fluid  respectively  [8].  This  shows  at  the  interface  between  the  magnet  and  fluid  there  is  a 
jump  discontinuity  in  the  magnetic  field  (assuming  that  the  magnetic  permeability  of  the  two 
materials  are  not  equal). 

The  model  presented  here  assumes  the  second  explanation  for  a  physical  interaction 
between  the  magnet  and  dipoles.  The  exact  force  caused  by  this  gradient  is  unknown,  but  it 
is  assumed  that  it  is  incredibly  short  ranged,  and  that  causes  a  force  of  attraction  such  that, 
when  multiplied  by  the  frictional  coefficient  between  the  particle  and  magnet,  leads  to  a 
frictional  force  that  is  substantially  larger  in  magnitude  as  compared  to  the  force  that  tends  to 
sweep  the  particle  along.  This  assumption  is  valid  since  the  force  tending  to  sweep  the 
particle  along  with  the  flow  is  a  function  the  flow  velocity  at  the  particle’s  location.  Since 
the  particle  is  small  (~5  microns)  and  resides  at  the  interface,  the  flow  velocity  is 
approximately  zero  (no  slip  condition).  In  other  word,  a  particle  that  happens  to  touch  the 
magnet  becomes  locked  in  place,  but  particles  in  the  stream,  away  from  the  wall,  do  not 
experience  this  force. 


16 


IV.  MODEL  FOR  INTERACTION 


A.  NEWTON’S  SECOND  LAW 


The  description  of  a  particle’s  motion  in  a  MR  fluid  can  be  detennined  using 
Newton’s  second  law  of  motion.  In  formulating  a  model  for  the  motion  of  an  arbitrary 
particle  i  apply  this  law  in  the  x  direction  as  follows 


d2x i 

=  m — d 
dt2 


(19) 


where  the  left  hand  side  of  the  equation  is  the  sum  of  all  of  the  forces  on  particle  i  in  the  x 
direction  and  m  is  the  particle’s  mass  (not  dipole  moment).  The  left  hand  side  includes 
the  dipole  interaction  force,  drag  force,  and  repulsive  forces  due  to  contact  with  other 
particles  and  the  walls.  Combining  equations  (7),  (10),  (13)  and  (16)  gives  the  following 
second  order  differential  equation  to  solve 


d2x. 
m — r1 

dt 2 


+  D 


(20) 


N 

where  D  =  6 jn]a  ,  and  Fix  =  FMix  +  ^  frepjj,x  +  fwu  ■ 

i*j 


Equation  (20)  can  be  solved  numerically,  in  its  current  form,  using  a  range  of 
techniques  (for  instance  a  Runge-Kutta  algorithm).  To  make  the  computations  more 
simple,  integrate  equation  (20)  over  a  sufficiently  small  time  step  r  such  that  the  term  on 
the  right  hand  side  can  be  assumed  constant.  This  gives  an  equation  for  the  change  in  the 
position  of  the  particle  in  the  x  direction  during  this  time  step  z  as  shown  below 


Ax  = 


,  Km 

D  D 


Dt 

l-e  m 


(21) 
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where  V0  is  the  velocity  of  the  particle  in  the  x  direction  at  the  beginning  of  the  time  step. 
This  equations  shows  that  if  r  is  chosen  such  that  it  is  several  orders  of  magnitude  less 
than  then  the  second  term  on  the  right  hand  side  can  be  ignored  and  therefore 


Ax  = 


FiS 

D 


and  similarly 


at 


D 


(22) 


(23) 


Az  = 


F.t 

iz 

D 


(24) 


for  the  y  and  z  components.  For  the  MR  fluids  analyzed  here  a  typical  value  of  yjj  is 

about  5E-7  sec'1  which  makes  the  time  step  on  the  order  of  IE-9  sec.  In  reality  this  will 
be  the  upper  limit  on  the  time  step.  Initially  a  much  smaller  time  step  will  be  used  in  the 
computer  simulation.  This  is  due  to  how  the  program  randomly  establishes  the  initial 
positions  of  the  particles  and  that  they  tend  to  overlap.  A  smaller  time  step  is  required  to 
“push”  the  particles  off  of  each  other  and  the  wall  without  destabilizing  the  computations 
with  excessively  large  positional  changes  at  each  time  step. 


B.  STATIC  FLUID  MODEL 


A  computational  algorithm  was  written  to  detennine  the  dynamic  motion  of  the 
particles  in  a  MR  fluid.  The  program  takes  user  inputs  for  the  length,  width  and  height  of 
the  MR  fluid  area,  the  number  of  particles  to  simulate,  the  magnitude  of  the  applied 
magnetic  field  (program  assumes  the  direction  is  in  the  negative  z  direction),  and  the 
number  of  time  steps  to  perform  the  algorithm.  The  program  then  randomly  distributes 
the  particles  inside  of  the  fluid  area.  Using  this  distribution  the  initial  spacing  between  all 
of  the  particles  is  computed  (the  values  ofxj,yj,zjj  and  n ).  Then  the  value  of  the  dipole 

strength  and  drag  coefficient  is  computed.  Using  the  spacing  between  particles  and  the 
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value  of  the  dipole  moments,  Fix ,  F.  ,  and  Fiz  are  calculated  for  every  particle.  If  the 

particle  is  near  the  upper  or  lower  wall  (a  distance  between  a  and  1.3a)  the  forces  are 
assumed  zero  for  the  reason  of  the  interaction  with  the  magnet/fluid  boundary  discussed 
in  Chapter  III.  Then  the  forces  are  computed  using  equations  (20).  A  time  step  is 

computed  based  on  the  value  of  ny^  as  described  in  the  previous  section.  The  updated 

position  of  every  particle  is  then  calculated  using  equations  (22-24)  and  this  new  position 
is  stored  and  plotted  graphically  if  desired.  This  process  is  repeated  for  every  time  step 
using  the  updated  positions  from  the  previous  time  step. 

As  discussed  above,  a  minor  issue  arises  in  the  initial  random  spacing,  especially 
at  higher  particle  densities.  Sometimes  the  particles  are  randomly  placed  such  that  two  or 
more  particles  are  spaced  where  the  distance  between  them  is  less  than  their  diameter 
length  apart  or  such  that  the  spacing  between  a  particle  and  a  wall  is  less  than  the 
particles  radius  length  apart.  This  is  not  physically  possible  since  the  particles  and  the 
wall  are  hard.  To  overcome  this,  the  time  step  chosen  for  the  first  10  time  steps  is  several 
orders  of  magnitude  less  than  what  is  used  for  the  remainder  of  the  computation.  This 
allows  the  repulsive  force  terms  to  “push”  the  particles  away  from  each  other  and  the 
wall  without  creating  an  abnormally  large  positional  change  that  would  eject  them  from 
the  MR  fluid  domain.  A  copy  of  the  computer  code  used  is  attached  in  Appendix  B  with 
a  more  detailed  discussion  as  to  the  inner  workings. 

C.  DYNAMIC  FLUID  MODEL 

The  programs  constructed  to  compute  the  dynamic  motion  of  particles  in  a 
dynamic  fluid  are  very  similar  to  the  one  for  the  static  case  with  a  few  alterations.  Two 
separate  programs  were  created,  one  for  pressure  driven  flow  and  the  other  for  shear 
driven  flow  (these  were  the  only  two  specific  flow  types  analyzed). 

In  the  pressure  driven  case  (parallel  flow  with  a  parabolic  velocity  distribution)  it 
is  assumed  that  the  flow  velocity  is  in  the  x  direction,  does  not  vary  in  the  x  and  y 
directions  and  varies  with  a  parabolic  distribution  in  the  z  direction.  The  user  inputs  the 
meanline  (maximum)  flow  and  the  program  computes  the  value  of  the  velocity  at  every 
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point  in  the  MR  fluid.  Using  this  velocity  distribution,  another  term  is  added  to  the  right 
hand  side  of  equation  (20)  to  account  for  the  drag  force  due  to  the  fluid  flow.  Equation 
(20)  now  becomes 


m 


d2x- 

~dd 


+  Ddxx  =  F+DU 

dt 


(25) 


where  Ui  is  the  flow  velocity  at  the  position  of  particle  i.  Note  that  since  the  flow  is  only 
in  the  x  direction  no  modification  is  required  for  the  equations  of  motion  in  the  y  and  z 
directions.  Applying  the  same  arguments  above  that  allowed  for  the  inertial  term  to  be 
ignored  allows  for  the  computation  of  the  change  in  the  position  of  the  particle  in  the 
same  manner  as  for  the  static  case. 

The  program  for  the  shear  driven  flow  (Couette  flow)  is  the  same  as  for  the 
pressure  driven  flow,  but  here  the  user  specifies  the  flow  velocity  at  the  upper  plate.  The 
program  then  calculates  the  velocity  at  all  other  points  in  the  fluid  and  simulates  the 
particle  motion  exactly  the  same  as  for  the  pressure  driven  flow. 
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V.  SIMULATION  RESULTS 


A.  STATIC  FLUID  SIMULATION 

The  first  qualitative  analysis  to  examine  is  the  effect  on  time  response  of  the  fluid 
as  a  function  of  particle  density.  The  results  shown  are  for  three  simulations  where  all 
parameters  are  held  constant  with  the  exception  of  particle  density.  The  various 
parameters  used  are  shown  in  the  below  table.  In  all  cases  the  size  of  the  rectangular 
volume  is  100  X  100  X  100  micrometers  with  hard  walls  bounding  the  area.  The 

magnetic  field  is  applied  in  the  negative  z  direction.  The  value  of  used  to  calculate 

the  time  step  has  a  value  of  1.744E-7  based  on  the  below  parameters. 


Number  of 
Particles 

Fluid  Viscosity 

Fluid 

Permeability 

Particle 

Permeability 

Applied 
Magnetic  Field 

Simulation  1 

40 

.25  Pa  s 

1.26E-6  N/A2 

.00377  N/A2 

200  kA/m 

Simulation  2 

70 

.25  Pa  s 

1.26E-6  N/A2 

.00377  N/A2 

200  kA/m 

Simulation  3 

too 

.25  Pa  s 

1.26E-6  N/A2 

.00377  N/A2 

200  kA/m 

Table  1.  Parameters  for  Simulations  1,  2  and  3 


The  simulations  were  conducted  out  for  100,000  time  steps  which  equates  to  a 
simulation  time  of  1.7  milliseconds.  The  figures  below  show  the  particle  microstructures 
at  various  times  for  the  above  simulations  in  both  a  3-D  and  top  down  view. 
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Figure  4.  Simulation  1,  Initial  Particle  Distribution,  3-D  View 
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Figure  5.  Simulation  1,  Initial  Particle  Distribution,  Top  Down  View 
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Figure  6.  Simulation  1,  Time  =  0. 16  milliseconds,  3-D  View 
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Figure  7.  Simulation  1,  Time  =  0. 16  milliseconds,  Top  Down  View 
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Figure  8.  Simulation  1,  Time  =  1.7  milliseconds,  3-D  View 
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Figure  9.  Simulation  1,  Time  =  1.7  milliseconds,  Top  Down  View 
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Figure  10.  Simulation  2,  Initial  Particle  Distribution,  3-D  View 


Figure  11.  Simulation  2,  Initial  Particle  Distribution,  Top  Down  View 
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Figure  12.  Simulation  2,  Time  =  0.16  milliseconds,  3-D  View 
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Figure  13.  Simulation  2,  Time  =  0.16  milliseconds,  Top  Down  View 
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Figure  14.  Simulation  2,  Time  =  0.86  milliseconds,  3-D  View 
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Figure  15.  Simulation  2,  Time  =  0.86  milliseconds,  Top  Down  View 
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Figure  16.  Simulation  2,  Time  =1.7  milliseconds,  3-D  View 


Figure  17.  Simulation  2,  Time  =1.7  milliseconds,  Top  Down  View 
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Figure  18.  Simulation  3,  Initial  Particle  Distribution,  3-D  View 


Figure  19.  Simulation  3,  Initial  Particle  Distribution,  Top  Down  View 
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Figure  20.  Simulation  3,  Time  =  0.16  milliseconds,  3-D  View 
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Figure  2 1 .  Simulation  3,  Time  =  0.16  milliseconds,  Top  Down  View 
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Figure  22.  Simulation  3,  Time  =  0.86  milliseconds,  3-D  View 
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Figure  23.  Simulation  3,  Time  =  0.86  milliseconds,  Top  Down  View 
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Figure  24.  Simulation  3,  Time  =1.7  milliseconds,  3-D  View 
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Figure  25.  Simulation  3,  Time  =1.7  milliseconds,  Top  Down  View 
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Qualitatively,  from  the  above  figures,  it  is  apparent  that  the  speed  at  which  the 
particles  form  chains  is  a  function  of  the  density  of  the  particles  in  the  fluid.  Examining 
the  two  extreme  cases  (Simulations  1  and  3),  much  longer  structures  with  fewer  smaller 
structures  are  evident  in  Figure  20  than  exist  in  Figure  6.  Therefore,  this  model 
demonstrates  the  experimentally  verified  fact  that  particle  density  is  a  factor  in  response 
time  of  the  MR  fluid  [3], 

The  second  qualitative  analysis  to  examine  is  the  effect  on  time  response  of  the 
fluid  as  a  function  of  applied  magnetic  field  strength.  The  results  shown  are  for  three 
simulations  where  all  parameters  are  held  constant  with  the  exception  of  the  magnetic 
field.  The  various  parameters  used  are  shown  in  Table  2.  In  all  cases  the  size  of  the 
rectangular  volume  is  100  X  100  X  100  micrometers  with  hard  walls  bounding  the  area. 

The  magnetic  field  is  applied  in  the  negative  z  direction.  The  value  of  'yjj  used  to 

calculate  the  time  step  has  a  value  of  1.744E-7  based  on  the  below  parameters. 


Number  of 
Particles 

Fluid 

Viscosity 

Fluid 

Permeability 

Particle 

Permeability 

Applied 

Magnetic 

Field 

Simulation  1 

70 

.25  Pa  s 

1.26E-6  N/A2 

.00377  N/A2 

150  kA/rn 

Simulation  2 

70 

.25  Pa  s 

1.26E-6  N/A2 

.00377  N/A2 

200  kA/rn 

Simulation  3 

70 

.25  Pa  s 

1.26E-6  N/A2 

.00377  N/A2 

250  kA/rn 

Table  2.  Parameters  for  Simulations  4,  5  and  6 
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Figure  26.  Simulation  4,  Initial  Particle  Distribution,  3-D  View 


Figure  27.  Simulation  4,  Initial  Particle  Distribution,  Top  Down  View 
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Figure  29.  Simulation  4,  Time  =  0.86  milliseconds,  Top  Down  View 
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Figure  30.  Simulation  4,  Time  =1.7  milliseconds,  3-D  View 
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Figure  3 1 .  Simulation  4,  Time  =1.7  milliseconds,  Top  Down  View 
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Figure  32.  Simulation  5,  Initial  Particle  Distribution,  3-D  View 
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Figure  33.  Simulation  5,  Initial  Particle  Distribution,  Top  Down  View 
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Figure  34.  Simulation  5,  Time  =  0.86  milliseconds,  3-D  View 
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Figure  35.  Simulation  5,  Time  =  0.86  milliseconds,  Top  Down  View 
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Figure  36.  Simulation  5,  Time  =1.7  milliseconds,  3-D  View 


Figure  37.  Simulation  5,  Time  =1.7  milliseconds,  3-D  View 
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Figure  38.  Simulation  6,  Initial  Particle  Distribution,  3-D  View 
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Figure  39.  Simulation  6,  Initial  Particle  Distribution,  Top  Down  View 
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Figure  40.  Simulation  6,  Time  =  0.86  milliseconds,  3-D  View 
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Figure  41.  Simulation  6,  Time  =  0.86  milliseconds,  Top  Down  View 
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Figure  42.  Simulation  6,  Time  =1.7  milliseconds,  3-D  view 
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Figure  43.  Simulation  6,  Time  =1.7  milliseconds,  Top  Down  View 
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Qualitatively,  from  the  above  figures,  it  is  apparent  that  the  speed  in  which  the 
particles  form  chains  is  a  function  of  the  strength  of  the  applied  magnetic  field. 
Examining  the  two  extreme  cases  (Simulations  4  and  6)  it  is  clear  that  the  structures  are 
much  closer  to  completion  in  Figure  40  than  those  in  Figure  28  (both  are  0.86 
milliseconds  into  the  simulation).  This  result  is  intuitive  based  on  the  fact  that  the 
strength  of  the  dipole  force  is  proportional  to  the  square  of  the  applied  magnetic  field 
strength  (Equations  1  and  3). 

B.  DYNAMIC  FLUID  SIMULATION 

The  dynamic  fluid  simulation  is  where  the  real  use  for  the  MR  model  is  realized. 
Most  MR  applications  in  industry  use  the  MR  effect  on  a  moving  fluid.  It  is  desirable  for 
the  dynamic  model  to  be  able  to  accurately  predict  the  microstructures  of  the  particles  in 
the  MR  fluid  and  from  there  predict  the  apparent  viscosity  and  shear  stress  of  the  MR 
fluid.  One  theory  for  predicting  the  shear  stress  has  already  been  developed  and  makes 
its  predictions  based  on  the  chain  density  in  the  MR  fluid  and  the  angles  the  chains  make 
in  relation  to  the  moving  fluid  [9]. 

A  simulation  of  a  MR  fluid  in  shear  was  run  in  order  to  show  that  these 
measurements  could  be  made  in  order  to  test  the  model  against  experimental  evidence 
and  to  use  the  model  to  help  design  MR  fluid  devices.  The  simulation  results  are  shown 
below  and  the  following  parameters  were  used. 


Number 

of 

Particles 

Fluid 

Viscosity 

Fluid 

Permeability 

Particle 

Permeability 

Applied 

Magnetic 

Field 

Velocity  of 
Top  Plate 

Simulation  7 

40 

.25  Pa  s 

1 .26E-6  N/A2 

.00377  N/A2 

250  kA/rn 

0.1  m/s 

Table  3.  Parameters  for  Simulation  7 


The  following  figures  show  the  distribution  of  particles  and  their  orientation  in  a 
dynamic  flow.  The  dimensions  for  the  area  were  changed  in  order  to  create  a  higher 
particle  density  without  adding  more  particles  to  the  volume.  In  this  case  the  size  of  the 
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rectangular  volume  is  50  X  100  X  50  micrometers  with  the  longer  direction  in  the 
direction  of  the  fluid  velocity.  The  particles  were  only  seeded  in  the  left  half  of  the 
volume  to  allow  the  fluid  to  push  the  particles  to  the  right  without  interference  from  the 

wall.  The  magnetic  field  is  applied  in  the  negative  z  direction.  The  value  of  nj/jj  used  to 

calculate  the  time  step  has  a  value  of  1.744E-7  as  was  the  case  for  the  above  simulations. 


x  10"5 


x  coordinate 

Figure  44.  Side  View  of  a  MR  Fluid  in  Shear 

In  the  above  figure  the  fluid  is  flowing  from  the  left  to  the  right  based  on  the  shear 
force  developed  due  to  the  top  plate  moving  to  the  left  at  0.1  m/s.  The  angle  that  the 
chains  make  with  relation  to  the  fluid  vary,  but  could  easily  be  measured  and  averaged. 
Examining  the  top  down  view  of  the  fluid  shown  below,  would  allow  for  the 
determination  of  the  chain  density. 
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Figure  45.  Top  View  of  the  MR  Fluid  in  Shear 

The  dynamic  model  allows  the  user  to  easily  input  common  variable  parameters 
in  a  program  in  order  to  detennine  the  dynamic  microstructure  in  two  common  flow 
conditions  (linear  flow  from  shear  and  parabolic  flow  from  a  pressure  gradient).  For 
other,  more  permanent  parameters  (particle  or  fluid  magnetic  penneability,  fluid 
viscosity,  etc.)  individual  lines  of  code  which  set  these  parameters  must  be  altered.  For  a 
more  detailed  description  of  what  parameters  the  user  inputs  when  running  the  code  and 
other  parameters  directly  set  in  the  code  see  Appendix  B. 

The  pressure  driven,  parabolic  velocity  profile  fluid  simulation  is  shown  below. 
Unfortunately  it  is  harder  to  determine  the  microstructure,  in  this  case  than  for  the  shear 
case.  The  red  line  inserted  on  the  below  figures  shows  one  chain  and  how  it  bulges  in  the 
middle  based  on  the  higher  flow  velocity  there.  The  parameters  for  this  simulation  are 
exactly  the  same  as  those  shown  in  Table  3,  except  the  flow  velocity  is  the  centerline 
flow,  not  the  flow  at  the  top  wall. 
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y  coordinate  z  coordinate 
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Figure  46.  Side  View  of  MR  Fluid  with  Parabolic  Flow 
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Figure  47.  Top  Down  View  of  MR  Fluid  with  Parabolic  Flow 
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VI.  CONCLUSION 


The  goal  of  the  model  developed  was  for  it  to  be  simple,  easy  to  use,  require  little 
empirical  data  (based  on  first  principles)  and  be  accurate.  The  model  satisfies  the 
requirement  to  be  simple.  It  uses  very  well  understood  laws  (Newton’s  Second  Law, 
dipole  interaction  force,  and  Stokes’  drag)  to  describe  the  forces  and  accelerations  of  the 
particles.  Through  a  mathematical  justification  it  ignores  the  inertial  mass  of  the  particles 
(being  dominated  by  the  viscous  forces)  to  simplify  the  calculations  even  further.  The 
physical  interaction  between  the  particles  and  the  walls  was  chosen  to  be  in  a  form  that 
would  balance  out  the  other  dominant  forces  in  a  way  that  was  short  ranged.  This 
technique  is  also  commonly  used  in  other  similar  types  of  models. 

The  programs  were  written  in  order  allow  multiple  simulations  with  a  minimum 
amount  of  work.  Instead  of  the  user  having  to  laboriously  input  every  parameter  required 
for  every  simulation  the  programs  only  require  the  user  to  input  a  small  number  of 
parameters  that  were  assumed  to  be  the  most  varied.  The  disadvantage  to  this  approach  is 
that  the  user  must  modify  the  code  in  order  to  change  parameters  such  as  particle  or  fluid 
magnetic  permeability,  fluid  viscosity,  or  particle  size.  However,  it  was  assumed  that 
these  parameters  would  not  be  changed  as  frequently  as  the  magnetic  field  strength,  size 
of  the  volume,  or  number  of  particles,  so  they  were  not  requested  by  the  program  for 
every  simulation. 

The  model  used  is  based  entirely  on  first  principles  so  no  empirical  data  is 
required  for  the  simulation  of  particles.  This  could  change  if  the  model  requires 
modification  in  order  to  accurately  predict  experimental  results.  For  instance,  in  some 
regimes,  the  Stokes’  approximation  may  no  longer  be  valid.  In  this  case  the  model 
should  be  modified  to  more  accurately  describe  the  shear  and  pressure  forces  on  a 
submerged  body  and  this  may  have  to  be  done  with  empirical  data.  It  is  the  hope  and 
belief  that  this  will  not  be  required,  but  it  is  a  possibility. 

As  far  as  accuracy  is  concerned  the  best  way  to  check  the  model  is  against 
carefully  controlled  laboratory  experiments.  Qualitatively  the  model  matches  observed 
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data  such  as  the  chain  formation  and  the  time  scale  in  which  these  structures  form  10. 
However,  there  are  several  assumptions  made  that  may  have  to  be  modified.  The  first 
assumption  that  limits  the  applicability  of  the  dynamic  model  is  the  characterization  of  the 
velocity  profile.  In  reality,  the  formation  of  particle  chains  alters  the  imposed  flow.  This  is 
how  MR  fluids  are  able  to  withstand  a  shear  and  alter  the  bulk  viscosity  of  the  fluid.  The 
model  presented  does  not  allow  for  the  changing  of  the  velocity  profile.  This  is  not  a  simple 
problem  to  solve,  because  the  only  way  to  determine  how  the  flow  changes  based  on  the 
particle  formations  and  how  the  particle  formations  vary  based  on  the  flow  is  to  numerically 
solve  to  sets  of  equations  simultaneously.  Either  Euler’s  equations  or  the  Navier-Stokes’ 
equations  must  be  solved  at  each  time  step  along  with  the  other  equations  to  determine  the 
forces  acting  on  the  particles.  This  would  require  the  integration  of  the  model  for  the 
dynamic  behavior  of  the  magnetic  particles  with  a  numerical  solver  for  fluid  dynamic.  The 
location  of  the  particles  at  each  time  step  would  represent  a  boundary  in  the  flow  that  would 
be  solved  with  a  flow  solver.  The  model  presented  here  does  not  attempt  to  perform  any  type 
of  alteration  of  the  flow  based  on  the  particle  dynamics.  Therefore,  the  model  is  not  good  for 
long  simulation  times,  but  it  is  still  valid  in  the  short  term  (before  the  initial  velocity  of  the 
fluid  is  altered).  In  other  words,  this  model  accurately  describes  the  initial  particle  dynamics, 
but  is  poor  in  the  limit  where  the  initial  flow  velocity  would  have  been  modified  by  the 
particle  structures. 

Another  area  requiring  more  detailed  study  is  in  the  interaction  between  the  particle 
chains  and  the  magnet  providing  the  magnetic  field.  As  discussed  earlier,  this  interaction  is 
well  understood  in  the  ER  case,  but  not  for  the  MR.  Models  for  ER  fluid  accurately  describe, 
based  on  the  interaction  with  the  chains  and  the  electrode,  how  particle  chains  will  merge  to 
form  even  larger  structures  around  one  second  after  the  field  is  applied.  These  larger 
structures  are  also  observed  in  MR  fluids,  but  the  method  in  which  they  form  is  being 
debated.  Since  it  is  doubted  that  the  chains  interact  with  the  magnet  in  the  same  way  that  the 
ER  fluid’s  chains  interact  with  the  electrodes,  the  same  process  is  not  occurring.  Some  have 
proposed  that  Brownian  motion  needs  to  be  included.  While  the  Brownian  motion  has  much 
smaller  interaction  energy  than  magnetism,  it  is  postulated  that  Brownian  interaction  could 
cause  the  chains  to  bulge  and  this  temporary,  minor  change  in  position  of  the  chain  could 
allow  for  nearby  chains  to  attract  this  chain.  Again,  this  is  a  postulation,  and  more  study  is 
required  for  the  particle-magnet  interaction. 
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APPENDIX  A.  DERIVATIONS  OF  EQUATION  4,  5,  AND  6 


Based  on  the  figure  below  define  fr  as  the  component  of  f..  in  the  radial 
direction  and  fg  as  the  component  of  ftj  in  the  angular  direction  as  shown. 


Figure  48.  Geometrical  Relationship  Between  Two  Particles 

After  dropping  the  subscripts  for  convenience  and  decomposing  fr  into  its 
Cartesian  components  gives 

fy  = sin(6')sin(^) 
f x  =  fr  sin(#)cos(^) 
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fz  =  fr  C0S(^) 


where  fy,  fx  and  fz  are  the  magnitudes  of  the  components  of  the  force  in  the  y,  x  and  z 
directions  respectively.  Similar  decompositions  of  fe  into  its  Cartesian  components 
gives 

fy  =fe  cos(0)sin(^) 
fx=  foCOS(0)COS(<f>) 
fz  =  -fe  sin(#) 

with  the  same  definitions  as  above.  Adding  the  components  together  gives 
fy  =  fr  sin(0)  sin(^)  +  fe  cos (0)  sin(^) 

fx  =  fr  sin(#)  C0S(^)  +  fe  C0S(#)  C0SW 

fz  =Xcos(6>)-/^sin(6>). 


Substituting  the  definition  of  fe  and  f0  from  equation  (3)  and  the  geometric  identities 


into  the  above  equation  obtains 
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Q  5z2  x 

f  =—  1 - 

Jx  4  2 

r  r  r 


f  =8-  3  5z~  £ 

/z  4  J  2 

r  r  r 


with  g  =  3m2///- 
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APPENDIX  B.  COMPUTER  CODE  AND  DESCRIPTION 


Below  is  the  computer  code  for  the  computation  of  the  static  flow  problem  and  a 
description  of  the  various  sections. 


%Thesis  program 
clear  all 

a=5*10A-6;  %m  radius  of  particle 
Vol=pi*aA3*4/3;  %Volume  of  particle 
mass=7 850*Vol ;  %kg  mass  of  particle 

tf=input (' Number  of  time  steps  '); 
height=input (' Height  of  volume  (micro  meter) '); 
height=height*10A-6;  %converts  to  meters 
length=input (' Length  of  volume  (micro  meter) '); 
length=length*10A-6;  %converts  to  meters 
width=input (' Width  of  volume  (micro  meter) '); 
width=width*10A-6; 

N=input (' Number  of  particles'); 

H=input (' Magnetic  Field  intensity  (kA/m) ');  %~200  kA/m 
H=H*1000; 

xinit=length*rand (N, 1 ) ;  %initial  x  dist  of  particles 
yinit=width*rand (N, 1 ) ;  %initial  y  dist  of  particles 
zinit=height*rand (N,  1 )  ;  %initial  z  dist  of  particles 

vis=.25;  %fluid  viscosity  [Pa*s] 
uf  =  1.257E-6;  %permeability  of  fluid 
up  =  .00377;  %permeability  of  particle 

m  =  ( 4 /3 ) *pi*H*aA3* (up-uf ) / (up+2 *uf ) ;  %magnetic  moment 
D=6*pi*vis*a;  %Stokes  drag  force  coefficient 
tcheck=mass/D;  %intrinsic  time  scale 
Q  =  3*mA2*uf ; 

ysys  =  zeros (tf,N);  %y  position  of  particles,  column  1  refers  to 
particle  1,  column  2  refers  to  particle  2,  ect 

xsys  =  zeros (tf,N);  %x  position  of  particles,  column  1  refers  to 
particle  1,  column  2  refers  to  particle  2,  ect 

zsys  =  zeros (tf,N);  %z  position  of  particles,  column  1  refers  to 
particle  1,  column  2  refers  to  particle  2,  ect 

delx  =  zeros (N,N);  %difference  in  the  x  position  between  the  particles 

dely  =  zeros (N,N);  %difference  in  the  y  position  between  the  particles 

delz  =  zeros (N,N);  %difference  in  the  z  position  between  the  particles 

r=zeros (N, N) ; 

Fmx=zeros (N, N) ; 

Fmy=zeros (N, N) ; 

Fpx=zeros (N,N) ; 

Fpy=zeros (N,N) ; 

Fpz=zeros (N,N) ; 
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Fmz=zeros (N,N)  ; 


ysys(l,:)  =  yinit ' ; 
xsys(l,:)  =  xinit'; 
zsys ( 1 , : )  =  zinit ' ; 


for  t  =  1 : tf 

x  =  xsys (t, : ) ; 
y  =  ysys (t,  : )  ; 
z  =  zsys (t, : ) ; 
for  i=l : N 

for  j  =1 : N 

delx  (i, j )  =  x (i) -x  ( j ) ; 
dely (i, j )  =  y (i) -y ( j ) ; 
delz  (i, j )  =  z  ( i ) - z  ( j ) ; 

r(i,j)  =  sqrt (delx (i, j ) A2+dely (i, j ) A2+delz  (i, j ) A2)  ; 
if  i== j 

Fmx (i, j ) =0; 

Fmy (i, j ) =0; 

Fpy (i, j ) =0 ; 

Fpx (i, j ) =0; 

Fpz (i, j ) =0; 

Fmz (i, j ) =0; 

else 


Fmx (i, j ) =-Q*delx (i,j)/r(i,j)A5*(5*(delz(i,j)/r(i,j))A2- 
1);  %force  on  particle  i  from  particle  j  in  x  direction  due  to  magnetic 
force 


Fpx (i, j ) =2*Q/ ( (2*a) A4) * (delx (i, j ) / (r(i,j) ) ) *exp (- 
12* ( (r (i,  j ) / (2*a) ) -1) ) ;  %force  on  particle  i  from  particle  j  in  x 
direction  due  to  physical  interaction  (collision) 

Fmy (i, j ) =-Q*dely (i, j)/r(i,j)A5*(5*(delz(i,j)/r(i,j))A2- 
1) ;  %force  on  particle  i  from  particle  j  in  y  direction  due  to  magnetic 
force 


Fpy (i, j ) =2*Q/ ( (2*a) A4) * (dely (i, j ) / (r(i,j) ) ) *exp (- 
12* ( (r (i, j ) / (2*a) ) -1) ) ;  %force  on  particle  i  from  particle  j  in  y 
direction  due  to  physical  interaction  (collision) 

Fmz  (i,  j  )  =-Q*delz  (i,j)  /r  (i,  j  )  A5*  (5*  (delz  (i,  j  )  /r  (i,  j  )  )  A2- 


3)  ; 


Fpz  ( i ,  j  )  =2  *Q/  (  ( 2  *  a )  A4)  *  (delz  (i,  j  )  /  (r  (i,  j  )  )  )  *exp  (- 
12*((r(i,j)/ ( 2  *  a ) ) - 1 ) ) ; 

Fx=Fmx+Fpx;  %total  force  in  x  direction 
Fy=Fmy+Fpy;  %total  force  in  y  direction 
Fz=Fmz+Fpz;  %total  force  in  z  direction 

end 


end 

Ftx=sum (Fx, 2 ) ;  %total  force  on  each  particle  in  x  direction 

Fty=sum (Fy, 2 ) ;  %total  force  on  each  particle  in  y  direction 

Ftz=sum (Fz , 2 ) ;  %total  force  on  each  particle  in  z  direction 

end 

for  q=l : N 

Ftz (q) =Ftz (q) +2*Q/ ( (2*a) A4) *exp (-30* (z (q)  / (2*a)  -  .5)  )  ; 

Ftz (q) =Ftz (q) -2*Q/ ( (2*a) A4) *exp (-30* ( (height-z (q) ) / (2*a) - 
.5));  %loop  incorporates  force  due  to  repulsion  of  wall 

Ftx (q) =Ftx (q) +2*Q/ ( (2*a) A4) *exp (-30* (x(q) / (2*a) -.5)  )  ; 
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Ftx (q) =Ftx (q) -2*Q/ ( (2*a) A4) *exp (-30* ( (length-x (q) ) / (2*a) - 

.  5 )  )  ; 

Fty (q) =Fty (q) +2*Q/((2*a)A4) *exp (-30*(y(q)/(2*a)-.5)); 

Fty (q) =Fty (q) -2*Q/ ( (2*a) A4) *exp (-30* ( ( width- y (q) ) / (2*a) - 

.  5 )  )  ; 

end 
if  t<10 

tau=t check/ 1E7 ; 
elseif  t<100 

tau=tcheck/ 1E4 ; 
elseif  t<1000 

tau=tcheck/ 1E2 ; 

else 

t au=t check/10; 

end 

xsys (t+1 , : ) =x+Ftx ' *tau/D; 
ysys (t+1 , : ) =y+Fty ' *tau/D; 
zsys (t+1 , : ) =z+Ftz ' *tau/D; 

plot3 (xsys (t+1,  :)  ,  ysys (t+1,  :) , zsys (t+1,  :) ,  'o',  'Markersize',25) ; 

axis  (  [0, length, 0, width, 0, height] ) 
pause ( . 0001 ) 
end 


The  first  section  clears  all  currently  stored  variables,  sets  particle  size  and 
calculates  the  mass  of  the  particles. 

%Thesis  program 
clear  all 

a=5*10A-6;  %m  radius  of  particle 
Vol=pi*aA3*4/3;  %Volume  of  particle 
mass=7850*Vol;  %kg  mass  of  particle 


The  next  section  is  for  the  user  to  input  the  number  of  time  steps,  volume 


dimensions,  number  of  particles,  and  magnetic  field  strength.  It  also  converts  these 


inputs  into  the  proper  units. 


tf=input (' Number  of  time  steps 
height=input (' Height  of  volume 
height=height*10A-6;  %converts 
length=input (' Length  of  volume 
length=length*10A-6;  %converts 
width=input (' Width  of  volume 
width=width*10A-6; 

N=input (' Number  of  particles') 
H=input (' Magnetic  Field  intens 


(micro  meter) ' ) ; 
to  meters 

(micro  meter) ' ) ; 
to  meters 
(micro  meter) ' ) ; 

ty  (kA/m) ');  %~200  kA/m 
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H=H*1000; 


The  next  section  of  code  randomly  decides  the  initial  positions  of  the 
particles  and  stores  them  in  the  variables  shown. 

xinit=length*rand (N, 1 ) ;  %initial  x  dist  of  particles 
yinit=width*rand (N, 1 ) ;  %initial  y  dist  of  particles 
zinit=height*rand (N, 1 ) ;  %initial  z  dist  of  particles 


The  next  section  set  some  parameters  of  the  MR  fluid  (magnetic 

permeability  of  the  particle  and  fluid  and  fluid  viscosity)  and  computes  the  magnetic 

dipole  moment,  Stokes’  drag  coefficient,  and  the  intrinsic  time  scale  (to  be  used  in 

determining  the  actual  time  scale  later  in  the  program). 

vis=.25;  %fluid  viscosity  [Pa*s] 
uf  =  1.257E-6;  %permeability  of  fluid 
up  =  .00377;  %permeability  of  particle 

m  =  ( 4 /3 ) *pi*H*aA3* (up-uf ) / (up+2 *uf )  ;  %magnetic  moment 
D=6*pi*vis*a;  %Stokes  drag  force  coefficient 
tcheck=mass/D;  %intrinsic  time  scale 
Q  =  3*mA2*uf ; 


The  next  section  preallocates  memory  for  the  matrices  that  will  be  used  to 
either  store  or  compute  the  motion  of  the  particles.  The  matrices  xsys,  ysys,  and  zsys  will 
store  the  actual  positions  of  the  particles.  The  columns  refer  to  the  individual  particles 
(the  first  column  is  the  position  of  particle  1 ,  the  second  column  refers  to  the  position  of 
particle  2,  etc.)  and  the  rows  refer  to  the  time  (the  first  row  is  the  initial  distribution,  the 
second  row  is  the  position  of  the  particles  after  one  time  step,  etc.).  The  matrices  delx, 
defy  and  delx  temporarily  store  the  differences  in  the  x,  y  and  z  direction  between 
particles.  The  index  of  the  matrix  determines  the  difference  in  position  of  which 
particles.  For  example,  delx(4,9)  is  the  difference  in  the  x  position  between  particles  4 
and  9.  These  matrices  get  written  over  after  each  time  step.  The  matrix  r  is  similar  to 
delx,  defy  and  delz  but  stores  the  difference  in  the  radial  direction  between  particles.  The 
matrices  Fmx,  Fmy,  and  Fmz  store  the  x,  y  and  z  components  of  the  force  between  two 
particles  due  to  their  magnetic  dipoles.  For  example,  Fmy(2,6)  is  the  dipole  force  in  the  y 
direction  between  particles  2  and  6.  Fpx,  Fpy,  and  Fpz  are  similar  to  Fmx,  Fmy,  and  Fmz 
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except  that  the  former  relate  to  the  physical  repulsive  force  due  to  the  particles  being  hard 
spheres. 

ysys  =  zeros (tf,N);  %y  position  of  particles,  column  1  refers  to 
particle  1,  column  2  refers  to  particle  2,  ect 

xsys  =  zeros (tf,N);  %x  position  of  particles,  column  1  refers  to 
particle  1,  column  2  refers  to  particle  2,  ect 

zsys  =  zeros (tf,N);  %z  position  of  particles,  column  1  refers  to 


particle  1,  column 

2  refers  to 

particle 

2 ,  ect 

delx  = 

zeros (N, N) ; 

%dif f erence 

in 

the 

X 

position 

between 

the 

particles 

dely  = 

zeros (N, N) ; 

%dif ference 

in 

the 

y 

position 

between 

the 

particles 

delz  = 

zeros (N, N) ; 

%dif ference 

in 

the 

z 

position 

between 

the 

particles 

r=zeros (N, N) ; 
Fmx=zeros (N, N) ; 
Fmy=zeros (N, N) ; 
Fmz=zeros (N,N) ; 
Fpx=zeros (N,N) ; 
Fpy=zeros (N,N) ; 
Fpz=zeros (N,N) ; 


The  next  section  takes  the  initial  particle  distribution  and  stores  this  data  in 
the  first  row  of  the  corresponding  position  matrices. 

ysys(l,:)  =  yinit'; 
xsys(l,:)  =  xinit'; 
zsys ( 1 , : )  =  zinit ' ; 


The  next  section  is  the  heart  of  the  program  that  calculates  the  dipole  and 
physical  interaction  between  particles.  These  nested  loops  calculate  first  the  delx,  defy, 
delz  and  r  matrices  for  the  current  time  step.  Then  set  the  diagonal  elements  of  the  force 
matrices  to  zero  (a  particle  does  not  interact  with  itself).  It  then  uses  the  equations 
discussed  in  the  thesis  to  calculate  all  of  the  other  elements  in  the  force  matrices. 

for  t  =  1 : tf 

x  =  xsys (t,  : )  ; 
y  =  ysys (t, : ) ; 
z  =  zsys (t, : ) ; 
for  i=l : N 

for  j  =1 : N 

delx  (i, j  )  =  x (i) -x  ( j  )  ; 
dely  (i, j  )  =  y ( i ) — y  ( j  )  ; 
delz  ( i , j  )  =  z (i) — z  (j)  ; 

r(i,j)  =  sqrt (delx (i, j ) A2+dely (i, j ) A2+delz (i, j ) A2) ; 
if  i== j 

Fmx (i, j ) =0; 

Fmy (i, j ) =0; 

Fpy (i, j ) =0; 
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Fpx (i, j ) =0; 

Fpz (i, j ) =0; 

Fmz (i, j ) =0; 

else 

Fmx (i, j ) =-Q*delx (i,j)/r(i,j)A5*(5*(delz(i,j)/r(i,j))A2- 
1) ;  %force  on  particle  i  from  particle  j  in  x  direction  due  to  magnetic 
force 

Fpx (i,j)=2*Q/((2*a)A4)*(delx(i,j)/(r(i,j))) *exp (- 
12* ( (r (i,  j ) / (2*a) ) -1) ) ;  %force  on  particle  i  from  particle  j  in  x 
direction  due  to  physical  interaction  (collision) 

Fmy (i, j ) =-Q*dely (i, j)/r(i,j)A5*(5*(delz(i,j)/r(i,j))A2- 
1);  %force  on  particle  i  from  particle  j  in  y  direction  due  to  magnetic 
force 


Fpy (i, j)=2*Q/((2*a)A4)*(dely(i,j)/(r(i,j))) *exp (- 
12* ( (r (i,  j ) / (2*a) ) -1) ) ;  %force  on  particle  i  from  particle  j  in  y 
direction  due  to  physical  interaction  (collision) 

Fmz (i, j ) =-Q*delz (i,j)/r(i,j)A5*(5*(delz(i,j)/r(i,j) 


3)  ; 


Fpz  (i,  j  )  =2*Q/  (  ( 2  *  a )  A4)  *  (delz  (i,  j  )  /  (r  (i,  j  )  )  )  *exp  (- 
12* ( (r (i, j ) / (2*a) ) -1) ) ; 

Fx=Fmx+Fpx;  %total  force  in  x  direction 
Fy=Fmy+Fpy;  %total  force  in  y  direction 
Fz=Fmz+Fpz;  %total  force  in  z  direction 

end 


A2- 


end 

Ftx=sum (Fx, 2 ) ; 
Fty=sum (Fy, 2 ) ; 
Ftz=sum (Fz , 2 ) ; 

end 


%total  force  on 
%total  force  on 
%total  force  on 


each  particle 
each  particle 
each  particle 


in  x  direction 
in  y  direction 
in  z  direction 


The  next  section  incorporates  a  loop  to  add  the  force  due  to  the  physical 
interaction  with  the  wall. 

for  q=l : N 

Ftz (q) =Ftz (q) +2*Q/ ( (2*a) A4) *exp (-30* (z (q)  / (2*a)  -  .5)  )  ; 

Ftz (q) =Ftz (q) -2*Q/ ( ( 2  *  a ) A4) *exp (-30* ( (height-z (q) ) /  (2*a)  —  .5) ) ; 
%loop  incorporates  force  due  to  repulsion  of  wall 

Ftx (q) =Ftx (q) +2*Q / ( ( 2  *  a ) A4) *exp (-30* (x(q) / (2*a) -.5) )  ; 

Ftx (q) =Ftx (q) -2*Q/ ( (2*a) A4) *exp (-30* ( (length-x (q))/(2*a)-.5)); 
Fty (q) =Fty (q)+2*Q/ ((2*a)A4) *exp (-30* (y (q) / (2*a) -.5) ) ; 

Fty (q) =Fty (q) -2*Q/ ( ( 2  *  a ) A4) *exp (-30* ( ( width- y (q) ) / (2*a) - . 5) )  ; 

end 


The  next  section  uses  the  intrinsic  time  scale  to  create  an  actual  time  scale 
based  on  the  time  that  the  program  is  simulating.  Extremely  small  time  scales  are  used 
initially  to  allow  the  initial  structures  to  begin  forming  and  then  the  time  scales  are 
enlarged  as  the  structures  become  closer  to  equilibrium. 
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if  t<10 


tau=t check/ 1E7 ; 
elseif  t<100 

tau=tcheck/ 1E4 ; 
elseif  t<1000 

tau=tcheck/ 1E2 ; 

else 

tau=t check/ 10; 

end 


The  final  section  calculates  the  new  position  of  every  particle,  stores  it  in 
the  appropriate  matrices  and  plots  the  positions  on  a  3-D  graph. 

xsys (t+1, : ) =x+Ftx ' *tau/D; 
ysys (t+1, : ) =y+Fty ' *tau/D; 
zsys (t+1 , : ) =z+Ftz ' *tau/D; 

plot3 (xsys (t+1,  :)  ,  ysys (t+1,  :) , zsys (t+1,  :) ,  'o',  'Markersize',25) ; 

axis  (  [ 0 , length, 0 , width, 0 , height] ) 
pause ( . 0001 ) 


Below  is  the  attached  code  for  the  shear  flow  problem.  This  code  is  very 
similar  to  the  static  code  above,  except  for  a  few  lines  described.  Sentences  in  red  are  the 
descriptions  of  the  changes  and  do  not  appear  in  the  code. 

%Thesis  program 
clear  all 

a=5*10A-6;  %m  radius  of  particle 
Vol=pi*aA3*4/3;  %Volume  of  particle 
mass=7 850*Vol ;  %kg  mass  of  particle 


tf=input (' Number  of  time  steps 
height=input (' Height  of  volume 
height=height*10A-6;  %converts 
length=input (' Length  of  volume 
length=length*10A-6;  %converts 
width=input (' Width  of  volume 
width=width*10A-6; 

U=input (' Velocity  of  top  plate 


(micro  meter) 
to  meters 

(micro  meter) 
to  meters 
(micro  meter) ' ) 

(m/s)');  I 


N=input (' Number  of  particles'); 
H=input (' Magnetic  Field  intensity 
H=H*1000; 


(kA/m)');  %~200  kA/m 


xinit=length*rand (N, 1 ) ;  %initial  x  dist  of  particles 
yinit=width*rand (N, 1 ) ;  %initial  y  dist  of  particles 
z ini t=height* rand (N, 1 ) ; 
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vis= . 25; 

uf  =  1.257E-6;  %permeability  of  fluid 
up  =  .00377;  %permeability  of  particle 

m  =  4*pi*H*aA3* (up-uf ) / (up+2*uf);  %magnetic  moment,  may  to  modify  H 
since  local  magnetic  field  may  not  be  applied  field  (see  paper 
micromechanical  model  for  MR  fluids) 

D=6*pi*vis*a;  %Stokes  drag  force  coefficient 
tcheck=mass/D;  %intrinsic  time  scale 
Q  =  3*mA2*uf ; 

ysys  =  zeros (tf,N);  %y  position  of  particles,  column  1  refers  to 
particle  1,  column  2  refers  to  particle  2,  ect 

xsys  =  zeros (tf,N);  %x  position  of  particles,  column  1  refers  to 

particle  1,  column  2  refers  to  particle  2,  ect 

zsys  =  zeros (tf,N); 

delx  =  zeros (N,N) ; 

dely  =  zeros (N,N); 

delz  =  zeros (N,N); 

r=zeros (N, N) ; 

Fmx=zeros (N, N) ; 

Fmy=zeros (N, N) ; 

Fpx=zeros (N,N) ; 

Fpy=zeros (N,N) ; 

Fpz=zeros (N,N) ; 

Fmz=zeros (N,N)  ; 

ysys(l,:)  =  yinit ' ; 
xsys(l,:)  =  xinit'; 
zsys ( 1 , : )  =  zinit ' ; 

for  t  =  1 : tf 

x  =  xsys (t,  : )  ; 
y  =  ysys (t, : ) ; 
z  =  zsys (t, : ) ; 
for  i=l : N 

for  j  =1 : N 

delx  (i, j  )  =  x (i) -x  ( j  )  ; 
dely  (i, j  )  =  y (i) -y  ( j  )  ; 
delz (i, j  )  =  z (i) — z  (j)  ; 

r(i,j)  =  sqrt (delx (i, j ) A2+dely (i, j ) A2+delz (i, j ) A2) ; 
if  i== j 

Fmx ( i , j ) =0 ; 

Fmy (i, j ) =0 ; 

Fpy (i, j ) =0; 

Fpx ( i , j ) =0 ; 

Fpz ( i , j ) =0 ; 

Fmz ( i , j ) =0 ; 

else 

Fmx (i, j ) =-Q*delx (i,j)/r(i,j)A5*(5*(delz(i,j)/r(i,j))A2- 
1) ;  %force  on  particle  i  from  particle  j  in  x  direction  due  to  magnetic 
force 
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Fpx (i,j)=Q/((2*a)A4)*(delx(i,j)/(r(i,j))) *exp (- 
12* ( (r (i, j ) / (2*a) ) -1) ) ;  %force  on  particle  i  from  particle  j  in  x 
direction  due  to  physical  interaction  (collision) 

Fmy (i, j ) =-Q*dely (i, j)/r(i,j)A5*(5*(delz(i,j)/r(i,j))A2- 
1) ;  %force  on  particle  i  from  particle  j  in  y  direction  due  to  magnetic 
force 


Fpy (i, j)=Q/((2*a)A4)*(dely(i,j)/(r(i,j))) *exp (- 
12* ( (r (i, j ) / (2*a) ) -1) ) ;  %force  on  particle  i  from  particle  j  in  y 
direction  due  to  physical  interaction  (collision) 

Fmz (i, j ) =-Q*delz (i,j)/r(i,j)A5*(5*(delz(i,j)/r(i,j) 


3)  ; 


Fpz  ( i ,  j  )  =2  *Q/  (  ( 2  *  a )  A4)  *  (delz  (i,  j  )  /  (r  (i,  j  )  )  )  *exp  (- 
12*((r(i,j)/ ( 2  *  a ) ) - 1 ) ) ; 

Fx=Fmx+Fpx;  %total  force  in  x  direction 
Fy=Fmy+Fpy;  %total  force  in  y  direction 
Fz=Fmz+Fpz ; 

end 


A2- 


end 

Ftx=sum (Fx, 2 ) ;  %total  force 
Fty=sum (Fy, 2 ) ;  %total  force 
Ftz=sum (Fz , 2 ) ;  %total  force 

end 

for  q=l : N 

if  z(q)<1.2*a  &&  z(q)>= 
Ftx (q) =  0; 

Fty (q) =0; 

Ftz (q) =  0; 

else 


on  each  particle 
on  each  particle 
on  each  particle 


a 


in  x  direction 
in  y  direction 
in  z  direction 


Ftz (q) =Ftz (q) +2*Q/ ( (2*a) A4) *exp (-30* (z (q) / (2*a) -.5) ) ; 

Ftz (q) =Ftz (q) -2*Q/ ( (2*a) A4) *exp (-30* ( (height-z (q) ) / (2*a) - 


.5));  %loop  incorporates  force  due  to  repulsion  of  wall 

Ftx (q) =Ftx (q) +2*Q/ ( (2*a) A4) *exp (-30* (x(q)  / (2*a)  -  .5)  )  ; 

Ftx (q) =Ftx (q) -2*Q/ ( (2*a) A4) *exp (-30* ( (length-x (q) ) / (2*a) - 


Fty (q) =Fty (q)+2*Q/ ( (2*a) A4) *exp (-30* (y (q) / (2*a) - . 5)  )  ; 

Fty (q) =Fty (q) -2*Q/ ( (2*a) A4) *exp (-30* ( ( width- y (q)  )  / (2*a)  - 


end 


end 


if  t<10 

tau=t check/ 1E7 ; 
elseif  t<100 

tau=tcheck/lE5 ; 
elseif  t<1000 

tau=t check/ 1E3 ; 

else 

tau=tcheck/100; 


end 
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xsys (t+1 , : ) =x+Ftx ' *tau/D; 
ysys (t+1, : ) =y+Fty ' *tau/D; 
zsys (t+1, : ) =z+Ftz ' *tau/D; 

plot3 (xsys (t+1,  :)  ,  ysys (t+1,  :)  , zsys (t+1,  :)  ,  'o',  'Markersize',25) 

axis  (  [ 0 , length, 0 , width, 0 , height] ) 
pause ( . 0001 ) 
end 
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